!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013
!
! This is for Y_k, R_k calculations
!
module jm_yk_mod;
use dbg_mod;
use jm_atom_mod;
use memory_mod;
use grid_mod;
use qdr_mod;
implicit none;
private;

! NOTE 12-Apr-2013: name any new public functions as yk_xxx()
public yk_lcr, yk_lcr_zk, yk_lcr_rk, yk_assert_lcr, test_jm_yk_mod;

!integer, parameter :: ATOM_ELEC_SPIN2 = 1; ! spin*2 for electron
real(DP), pointer, save :: WORK_VEC(:) => NULL();

contains;!=====================================================================
!___________________________________________________________
subroutine test_jm_yk_mod();
  character(LEN=64) :: name = "test_jm_yk_mod";
  call dbg_in(name);   
  call test_yk_assert_lcr_N221();
  call test_yk_assert_lcr_N1001();  
  call dbg_out(name); 
end subroutine
!___________________________________________________________
SUBROUTINE test_yk_assert_lcr_N221(); 
  type(GRID_LCR_T) :: lcr;
  character(LEN=64) :: name = "test_yk_assert_lcr_N221";
  call dbg_in(name);
  call grid_lcr_init(lcr);
  lcr%xmin = -D4;  lcr%rmin = D0;  lcr%nx = 221;
  lcr%xstep = D1/D16;
  call grid_lcr_xstep(lcr);
  call qdr_lcr(lcr);
  call yk_assert_lcr(lcr, EPS12, D5 * EPS8, EPS9);
end SUBROUTINE  
!___________________________________________________________
SUBROUTINE test_yk_assert_lcr_N1001(); 
  type(GRID_LCR_T) :: lcr;
  character(LEN=64) :: name = "test_yk_assert_lcr_N1001";
  call dbg_in(name);
  call grid_lcr_init(lcr);
  lcr%xmin = -D5;  lcr%rmax = D100;  lcr%nx = 1001;
  call grid_lcr(lcr);
  call qdr_lcr(lcr);
  call yk_assert_lcr(lcr, EPS13, EPS12, EPS13);
end SUBROUTINE 
!___________________________________________________________
SUBROUTINE yk_lcr(res, wf, wf2, k, lcr)
  real(DP), intent(OUT) :: res(:); ! result vector
  real(DP), intent(IN) :: wf(:), wf2(:);
  integer, intent(IN) :: k;
  type(GRID_LCR_T), intent(IN) :: lcr;
  character(LEN=64) :: name = "yk_lcr";
  call dbg_in(name);
  !call manage_work_zk(wf, wf2, lcr);  
  call calcZkV2(res, wf, wf2, k, lcr); 
  call calcYkV2(res, wf, wf2, k, lcr);  
end SUBROUTINE
!___________________________________________________________
!// from F.Fisher "The hartree-fock method for atoms", p18
!//                     oo oo
!// R^k(a, b; a2, b2) = I  I  dr ds    P(a; r) P(a2; r) U^k(r,s) P(b; s) P(b2; s)
!//                     0  0
!// = INTL dr P(a; r) P(a2; r) 1/r Y(b, b2, K, r); NOTE 1/r!!!
!//  U^k(r,s) = y^k/x^(k+1)   y = min(r,s)   x = max(r,s)
real(DP) function yk_lcr_rk(a, b, a2, b2, k, lcr) result(res);
  real(DP), intent(IN) :: a(:), b(:), a2(:), b2(:);
  integer, intent(IN) :: k;
  type(GRID_LCR_T), intent(IN) :: lcr;
  character(LEN=64) :: name = "yk_lcr_rk";
  call dbg_in(name);
  call manage_work_vec(b, b2, lcr);    
  call yk_lcr(WORK_VEC, b, b2, k, lcr);
  res = qdr_intgrl_lcr2DivR(lcr, WORK_VEC, a, a2);   
end function
!___________________________________________________________
SUBROUTINE manage_work_vec(wf, wf2, lcr)
  real(DP), intent(IN) :: wf(:), wf2(:);
  type(GRID_LCR_T), intent(IN) :: lcr;
  integer :: mini, maxi, minw, maxw, n, nw;
  if (ASSOCIATED(WORK_VEC)) then
    call vec_range(mini, maxi, wf, n);
    call vec_range(minw, maxw, WORK_VEC, nw);
    if (nw >= n) return !enough space
    call mem_del_vec(WORK_VEC);
  endif
  WORK_VEC => mem_vec_from(wf);
end SUBROUTINE
!___________________________________________________________
SUBROUTINE yk_lcr_zk(res, wf, wf2, k, lcr)
  real(DP), intent(OUT) :: res(:); ! result vector
  real(DP), intent(IN) :: wf(:), wf2(:);
  type(GRID_LCR_T), intent(IN) :: lcr;
  integer, intent(IN) :: k;
  call calcZkV2(res, wf, wf2, k, lcr); 
end SUBROUTINE
!___________________________________________________________
SUBROUTINE calcYkV2(res, wf, wf2, k, lcr)
  real(DP), intent(INOUT) :: res(:); ! in zk, out:result yk
  real(DP), intent(IN) :: wf(:), wf2(:);
  integer, intent(IN) :: k;
  type(GRID_LCR_T), intent(IN) :: lcr;
  integer :: i, mini, maxi, nr, K1;
  real(DP) :: t90, t114, t34, C, f1, f2, f3, f4, f5, ai, a1, a2, a3, rM1;
  character(LEN=64) :: name = "calcYkV2";
  call dbg_in(name);
  
  call vec_range(mini, maxi, lcr%rv, nr);
  call dbg_assert_zero("lcr%rv(mini)==0", lcr%rv(mini));

  K1 = k + 1;
  C = real(2 * k + 1, DP);  
  !! the following t90, t114, t34 are diff from lcr%t90... when k.not.0
  t90 = lcr%xstep / D90 * C;!!!!// NOTE!!! *C
  t114 = D114 * t90;
  t34 = D34 * t90;

!  int last = res.length - 1;  // starting from last
  f1 = lcr%crDivR(maxi)   * res(maxi);
  f2 = lcr%crDivR(maxi-1) * res(maxi-1);
  f3 = lcr%crDivR(maxi-2) * res(maxi-2);
  f4 = lcr%crDivR(maxi-3) * res(maxi-3);
! for (int i = last - 2; i >= 2; i--) {
  do i = maxi - 2, 3, -1
!   f5 = crDivR[i - 2] * res[i - 2];
    f5 = lcr%crDivR(i - 2) * res(i - 2);
!   rM1 = r[i - 1]; !// NOTE K1 not K below
    rM1 = lcr%rv(i - 1); !// NOTE K1 not K below
!   a3 = Mathx.pow(rM1 * divR[i + 2], K1) * h90;
    a3 = (rM1 * lcr%divR(i + 2))**K1 * t90;
!   a2 = Mathx.pow(rM1 * divR[i + 1], K1);
    a2 = (rM1 * lcr%divR(i + 1))**K1;
!   a1 = Mathx.pow(rM1 * divR[i + 0], K1) * h114;
    a1 = (rM1 * lcr%divR(i + 0))**K1 * t114;
!   ai = Mathx.pow(rM1 * divR[i - 2], K1) * h90;
    ai = (rM1 * lcr%divR(i - 2))**K1 * t90;
!//      YK(M) = YK(M+2)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)        AATK4078
!   res[i - 1] = res[i + 1] * a2 + (a1 * f3 + a34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
    res(i - 1) = res(i + 1) * a2 + (a1 * f3 + t34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
    f1 = f2;
    f2 = f3;
    f3 = f4;
    f4 = f5;
  enddo
  res(mini) = D0;
end SUBROUTINE
!___________________________________________________________
SUBROUTINE calcZkV2(res, wf, wf2, k, lcr)
  real(DP), intent(OUT) :: res(:); ! result vector
  real(DP), intent(IN) :: wf(:), wf2(:);
  integer, intent(IN) :: k;
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP) :: f1, f2, f3, f4, f5, a, b, intgl, norm, corr; 
  real(DP) :: r2, r3, rk2, rk3, f21, f31, hb, rkn2, rkn3, ovR; 
  real(DP) :: ai, a1, a2, a3; 
  real(DP), pointer :: powR(:);
  integer :: mini, maxi, n, i;
  character(LEN=64) :: name = "calcZkV2";
  call dbg_in(name);

  !call vec_dbg("lcr%cr2=", lcr%cr2);
  !call vec_dbg("wf=", wf);
  !call vec_dbg("wf2=", wf2);
  call vec_range(mini, maxi, res, n);
  res(mini) = D0;
  f1 = lcr%cr2(mini) * wf(mini) * wf2(mini);        !call dbg_dp("f1=", f1);
  call dbg_assert_zero("f1==0", f1);
  f2 = lcr%cr2(mini+1) * wf(mini+1) * wf2(mini+1);  !call dbg_dp("f2=", f2);
  f3 = lcr%cr2(mini+2) * wf(mini+2) * wf2(mini+2);  !call dbg_dp("f3=", f3);
  f4 = lcr%cr2(mini+3) * wf(mini+3) * wf2(mini+3);  !call dbg_dp("f4=", f4);

!  // NOTE!!! Integrating by r (not X)
!  // y(r) = a r^b + c
!  // f1 = c
!  // f2 = a r2^b + c
!  // f3 = a r3^b + c
!  // f2 - f1 = a r1^b
!  // f3 - f1 = a r2^b
!  // (f3-f1)/(f2-f1)= (r3/r2)^b;
  r2 = lcr%rv(mini+1);
  r3 = lcr%rv(mini+2);
  rk2 = r2**k;
  rk3 = r3**k;
  f21 = f2 * rk2 / lcr%cr(mini+1);
  f31 = f3 * rk3 / lcr%cr(mini+2);
  b = LOG(f31 / f21) / LOG(r3 / r2);    !call dbg_dp("b=", b);
  call dbg_assert_true("b > 0", b > 0);
  call dbg_assert_not_zero("b.not.==0", b);
  hb = r2**b;
  a = f21 / hb;
  intgl = a *hb / (b+D1) * r2;   !call dbg_dp("res[0]=", res(1));
  res(mini+1) = intgl / rk2;     !call dbg_dp("res[1]=", res(2));
!
   hb = r3**b;
  intgl = a *hb / (b+D1) * r3;
  res(mini+2) = intgl / rk3;     !call dbg_dp("res[2]=", res(3));
!
!  // START DEBUG CODE ====================================
!  //
!  // Simpson rule
!  // y(x) = a x^2 + b x + c
!  // f1 = c
!  // f2 = a h^2 + b h + c
!  // f3 = 4 a h^2 + 2 b h + c
!  //  double a = (f3 - 2. * f2 - c) / (2. * h * h);
!  //  double b = (4. * f2 - f3 + 3c) / (2. * h);
!  // int_0^h y(x) dx = a / 3. * h^3 + b h^2 / 2. + c h;
!  //  a = (f3 - 2. * f2 + f1) / 2.;
!  //  b = (4. * f2 - f3 - 3. * f1) / 2.;
!  //  intgl = (a / 3.  + b / 2. + f1) * h;
!  //  res[1] = res[0] + intgl;
!  //  res[1] = h2 * f2;   // trapezoidal
!  //
!  //  // Simpson rule
!  //  a = Mathx.pow(r[1] * divR[2], k);
!  //  //YK(3) = YK(1)*A2 + H3*(F3 + D4*A*F2 + A2*F1)
!  //  res[2] = h3*(f3 + 4.*a*f2);// NOTE!!! r[0]=Zk[0]=a2=f1=0 in the LCR grid
!  //
!  // END DEBUG CODE ====================================

!  for (int i = 4; i < n; i++) {
  do i = 4, n-1  
!   f5 = cr2[i] * wf.get(i) * wf2.get(i);//      F5 = (RR(M)*P(M,I))*P(M,J)                                        AATK4125
    f5 = lcr%cr2(mini+i) * wf(mini+i) * wf2(mini+i); !//      F5 = (RR(M)*P(M,I))*P(M,J)                                        AATK4125

    !// Zim = (ri/rim)^k Zi + I^im_i FF y^2 * (r/rim)^k dx;  where ri=r_i, rim=r_(i+m)
    !// the same as in Froese-Fischer but replacing
    !//  exp(-mhk) with (ri/rim)^k;   and
    !//  r*exp(k(x-xim)) with  y^2*(r/rim)^k

    !// Replace exp(-mhk) with (ri/rim)^k;  where m=2
    !//      EH = DEXP(-h)   // from SUBROUTINE INIT
    !//      A = EH**K
!   ovR = divR[i - 1];
    ovR = lcr%divR(mini+i-1);
    a3 = (lcr%rv(mini+i-4)*ovR)**k * lcr%h90; !//      A3 = A2*A*H90                                                     AATK4112
    a2 = (lcr%rv(mini+i-3)*ovR)**k;       !//      A2 = A*A                                                          AATK4110
    a1 = (lcr%rv(mini+i-2)*ovR)**k * lcr%h114;!//      AN = 114.D0*A*H90                                                 AATK4114
    ai = (lcr%rv(mini+i-0)*ovR)**k * lcr%h90; !//      AI = H90/A                                                        AATK4113

    !// From Froese Fischer etal Computational atomic structure 2000
    !// DELTA^2 F = F2 - 2*F3 + F4
    !// DELTA^4 F = F1 - 4*F2 + 6*F3 - 4*F4 + F5
    !// INTL F(x)dx = h(F3 + DELTA^2 F/3 - DELTA^4 F/90) + O(h^7) =
    !//  = h/90 [114*F3 + 34*(F2+F4) - F1 - F5]
    !//
    !//      YK(M-1) = YK(M-3)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)
    !res[i-1] = res[i-3] * a2 + (a1 * f3 + a34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
    res(mini+i-1) = res(mini+i-3) * a2 + (a1 * f3 + lcr%a34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
    !// NOTE: missing i=n-1
    f1 = f2;
    f2 = f3;
    f3 = f4;
    f4 = f5;
  enddo
  !// NOTE: missing i=n-1
  norm = D0;
  rkn2 = D1;
  rkn3 = D1;
  if (k == 0) then  !// correct from last two points
!   norm = quadr.calcInt(wf, wf2);
    norm = qdr_intgrl_lcr2(lcr, wf, wf2);
  else 
!   norm = quadr.calcInt(wf, wf2, quadr.getPowR(k));
    powR => grid_lcr_get_pow_r(lcr, k);
    norm = qdr_intgrl_lcr3(lcr, wf, wf2, powR);
!   rkn2 = 1. / quadr.getPowR(k).get(n - 2);
    rkn2 = D1 / powR(mini+n-2);
!   rkn3 = 1. / quadr.getPowR(k).get(n - 3);
    rkn3 = D1 / powR(mini+n-3);
  endif
  corr = norm * rkn2 - res(mini+n-2);      !//log.info("corr=", corr);
! for (int i = n-2; i > 0; i -= 2) { // do not correct i=0
  do i = n-2, 1, -2 ! // do not correct i=0
    res(mini+i) = res(mini+i) + corr;
  enddo
  corr = norm * rkn3 - res(mini+n-3);!             //log.info("corr=", corr);
! for (int i = n-3; i > 0; i -= 2) { // do not correct i=0
  do i = n-3, 1, -2  ! // do not correct i=0
    res(mini+i) = res(mini+i) + corr;
  enddo
  !//last point
  a1 = (lcr%rv(mini+n-2) * lcr%divR(mini+n-1))**k;
  res(mini+n-1) = a1 * res(mini+n-2);
end SUBROUTINE
!___________________________________________________________
!
!! TESTS ===================================================
!___________________________________________________________
! call this to test an lcr-grid
! Rk values are from p236 of C.F.Fischer
!public void testRkLcr() throws Exception {
!    double FIRST = -4;
!    int NUM_STEPS = 220;
!    double STEP = 1. / 16.;
!    StepGrid x = new StepGrid(FIRST, NUM_STEPS, STEP);
!    TransLcrToR xToR = new TransLcrToR(x);
!    Vec r = xToR;
!//      QuadrPts5 w = new QuadrPts5(x);
!    WFQuadrLcr wCR = new WFQuadrLcr(x);
!
!    FuncVec f = WfFactory.makeP1s(r, 1.);
!    f.multSelf(xToR.getDivSqrtCR());
!    double res = 0, rMm = 0;
!    res = wCR.getWithCR2().calc(f, f);
!    assertEquals(0, Math.abs(res - 1), 6e-13);
!
!    FuncVec f2 = WfFactory.makeP1s(r, 1.);
!    f2.multSelf(xToR.getDivSqrtCR());
!    f2.mult(2.); // NOTE! Checking diff norm
!    res = wCR.getWithCR2().calc(f2, f2);
!    assertEquals(0, Math.abs(res - 4), 3e-12);
!
!    FuncVec T = WfFactory.makeY_0_1s(r); // valid
!    FuncVec Y = new YkLcr(wCR, f, f, 0).calcYk();
!//    assertEquals(0, Math.abs(T.distSLOW(Y)), 2e-9);
!    assertEquals(0, Math.abs(DistMaxAbsErr.distSLOW(T, Y)), 2e-9);
!    res = RkLcr.calc(wCR, f, f, f, f, 0);
!    // from p236 of C.F.Fischer
!    assertEquals(0, Math.abs(res - 5. / 8.), 5e-11);  //f95
!
!    res = RkLcr.calc(wCR, f2, f, f, f, 0);  // note f2
!    assertEquals(0, Math.abs(res - 2. * 5. / 8.), 9e-11);
!
!    res = RkLcr.calc(wCR, f, f2, f, f, 0);  // note f2
!    assertEquals(0, Math.abs(res - 2. * 5. / 8.), 9e-11);
!
!    res = RkLcr.calc(wCR, f, f, f2, f, 0);  // note f2
!    assertEquals(0, Math.abs(res - 2. * 5. / 8.), 9e-11);
!
!    res = RkLcr.calc(wCR, f, f, f, f2, 0);  // note f2
!    assertEquals(0, Math.abs(res - 2. * 5. / 8.), 9e-11);
!
!    HkMm hkMm  = new HkMm(wCR, f, f, f, f, 0);
!    rMm = hkMm.calcPot2();
!    assertEquals(0, Math.abs(rMm - 5. / 8.), 2e-10); // todo: was 5e-11???????
!
!    FuncVec f2s = WfFactory.makeP2s(r, 1.);
!    f2s.multSelf(xToR.getDivSqrtCR());
!    res = wCR.getWithCR2().calc(f2s, f2s);
!    assertEquals(0, Math.abs(res - 1), 2e-13);
!
!    res = RkLcr.calc(wCR, f, f2s, f, f2s, 0);
!    assertEquals(0, Math.abs(res - 17. / 81), 4e-11);
!
!    res = RkLcr.calc(wCR, f2s, f, f2s, f, 0);
!    assertEquals(0, Math.abs(res - 17. / 81), 4e-11);
!    res = RkLcr.calc(wCR, f, f2s, f2s, f, 0);
!    assertEquals(0, Math.abs(res - 16. / 729), 5e-11);
!    res = RkLcr.calc(wCR, f2s, f, f, f2s, 0);
!    assertEquals(0, Math.abs(res - 16. / 729), 5e-11);
!
!    res = RkLcr.calc(wCR, f2s, f2s, f2s, f2s, 0);
!    hkMm  = new HkMm(wCR, f2s, f2s, f2s, f2s, 0);
!    rMm = hkMm.calcPot2();
!    assertEquals(0, Math.abs(res - 77. / 512), 4e-10);
!    assertEquals(0, Math.abs(rMm - 77. / 512), 3e-10);
!
!    FuncVec f2p = WfFactory.makeP2p(r, 1.);
!    f2p.multSelf(xToR.getDivSqrtCR());
!    res = wCR.getWithCR2().calc(f2p, f2p);
!    assertEquals(0, Math.abs(res - 1), 7e-14);
!
!    res = RkLcr.calc(wCR, f2p, f, f2p, f, 0);
!    assertEquals(0, Math.abs(res - 59. / 243), 3e-11);
!    res = RkLcr.calc(wCR, f, f2p, f, f2p, 0);
!    assertEquals(0, Math.abs(res - 59. / 243), 3e-11);
!    res = RkLcr.calc(wCR, f2p, f, f, f2p, 1);
!    assertEquals(0, Math.abs(res - 112. / 2187), 2e-10);
!    res = RkLcr.calc(wCR, f2s, f2p, f2s, f2p, 0);
!    assertEquals(0, Math.abs(res - 83. / 512), 2e-10);
!    res = RkLcr.calc(wCR, f2p, f2s, f2p, f2s, 0);
!    assertEquals(0, Math.abs(res - 83. / 512), 2e-10);
!    res = RkLcr.calc(wCR, f2s, f2p, f2p, f2s, 1);
!    assertEquals(0, Math.abs(res - 45. / 512), 9e-10);
!    res = RkLcr.calc(wCR, f2p, f2p, f2p, f2p, 0);
!    assertEquals(0, Math.abs(res - 93. / 512), 1e-10);
!    res = RkLcr.calc(wCR, f2p, f2p, f2p, f2p, 1);
!    assertEquals(0, Math.abs(res - 0.12044270784428178), 1e-16); // at x=-4, 1/16, 220
!    res = RkLcr.calc(wCR, f2p, f2p, f2p, f2p, 2);
!    assertEquals(0, Math.abs(res - 0.08789062488574255), 1e-16); // at x=-4, 1/16, 220
!  }
SUBROUTINE yk_assert_lcr(lcr, normErr, yzErr, rkErr); 
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: normErr, yzErr, rkErr;
  real(DP), pointer :: f(:), f2(:), f2p(:), tk(:), ck(:); 
  real(DP) :: res, rk;
  integer :: n, L;
  character(LEN=64) :: help;
  character(LEN=64) :: name = "yk_assert_lcr";
  call dbg_in(name);

  f => mem_vec_from(lcr%rv);
  f2 => mem_vec_from(lcr%rv);
  f2p => mem_vec_from(lcr%rv);
   
!  // 1s-1s tests
  call atom_wfv_h1s(f, lcr%rv, D1);    call vec_dbg("h1s=", f);
  call grid_lcr_wf(f, lcr);            call vec_dbg("h1s_lcr=", f);
  res = qdr_intgrl_lcr2(lcr, f, f);    call dbg_dp("res=", res);
  call dbg_assert_abs_LE("< P1s | P1s > = ", res - D1, normErr);
   
  tk => mem_vec_from(lcr%rv);
  call atom_z1v_1s(tk, lcr%rv);        call vec_dbg("tZk=", tk);
  ck => mem_vec_from(lcr%rv);
  call yk_lcr_zk(ck, f, f, 1, lcr);    call vec_dbg("cZk=", ck);
  call vec_assert_abs_diff_LE("Z_k1_1s = ", tk, ck, yzErr); 

  help = "Z_k0_1s = ";
  call atom_z0v_1s(tk, lcr%rv);        call vec_dbg("tZk=", tk);
  call yk_lcr_zk(ck, f, f, 0, lcr);    call vec_dbg("cZk=", ck);
  call vec_assert_abs_diff_LE(TRIM(help), tk, ck, yzErr); 

  help = "Y_k0_1s = ";
  call atom_y0v_1s(tk, lcr%rv);        call vec_dbg("tYk=", tk);
  call yk_lcr(ck, f, f, 0, lcr);       call vec_dbg("cYk=", ck);
  call vec_assert_abs_diff_LE(TRIM(help), tk, ck, yzErr); 

  help = "R_k0_1s = ";
  res = yk_lcr_rk(f, f, f, f, 0, lcr); call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - D5/D8, rkErr);  

!  // 1s-2s tests
  n = 2; L = 0;
  call atom_wfv_hnL(f2, lcr%rv, n, L, D1);  call vec_dbg("h2s=", f2);
  call grid_lcr_wf(f2, lcr); call vec_dbg("h2s_lcr=", f2);
  res = qdr_intgrl_lcr2(lcr, f2, f2);  call dbg_dp("res=", res);
  call dbg_assert_abs_LE("< P2s | P2s > = ", res - D1, normErr);
  res = qdr_intgrl_lcr2(lcr, f, f2);  call dbg_dp("res=", res);
  call dbg_assert_abs_LE("< P1s | P2s > = ", res, normErr);

  call atom_y0v_2s(tk, lcr%rv);  call vec_dbg("tYk=", tk);
  call yk_lcr(ck, f2, f2, 0, lcr);  call vec_dbg("cYk=", ck);
  call vec_assert_abs_diff_LE("Y_k0_2s = ", tk, ck, yzErr); 

  help = "R_k0_1s_2s_1s_2s = ";  rk = D17 / 81._DP;
  res = yk_lcr_rk(f, f2, f, f2, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2, f, f2, f, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k0_1s_2s_2s_1s = ";  rk = D16 / 729._DP;
  res = yk_lcr_rk(f, f2, f2, f, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2, f, f, f2, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k0_2s= ";  rk = 77._DP / 512._DP;
  res = yk_lcr_rk(f2, f2, f2, f2, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

!  // 2p, k=0
  n = 2; L = 1;
  call atom_wfv_hnL(f2p, lcr%rv, n, L, D1);    call vec_dbg("h2p=", f2p);
  call grid_lcr_wf(f2p, lcr);                  call vec_dbg("h2p_lcr=", f2p);
  res = qdr_intgrl_lcr2(lcr, f2p, f2p);        call dbg_dp("res=", res);
  call dbg_assert_abs_LE("< P2p | P2p > = ", res - D1, normErr);

  call atom_y0v_2p(tk, lcr%rv);                call vec_dbg("tY_k0_2p=", tk);
  call yk_lcr_zk(ck, f2p, f2p, 0, lcr);        call vec_dbg("cZ_k0_2p=", ck);
  call yk_lcr(ck, f2p, f2p, 0, lcr);           call vec_dbg("cY_k0_2p=", ck);
  call vec_assert_abs_diff_LE("Y_k0_2p = ", tk, ck, yzErr); 

  help = "R_k0_1s_2p_1s_2p= ";  rk = 59._DP / 243._DP;
  res = yk_lcr_rk(f, f2p, f, f2p, 0, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2p, f, f2p, f, 0, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k1_1s_2p_2p_1s= ";  rk = 112._DP / 2187._DP;
  res = yk_lcr_rk(f, f2p, f2p, f, 1, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2p, f, f, f2p, 1, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k0_2s_2p_2s_2p= ";  rk = 83._DP / 512._DP;
  res = yk_lcr_rk(f2, f2p, f2, f2p, 0, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2p, f2, f2p, f2, 0, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k1_2s_2p_2p_2s= ";  rk = 45._DP / 512._DP;
  res = yk_lcr_rk(f2, f2p, f2p, f2, 1, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);
  res = yk_lcr_rk(f2p, f2, f2, f2p, 1, lcr);     call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k0_2p= ";  rk = 93._DP / 512._DP;
  res = yk_lcr_rk(f2p, f2p, f2p, f2p, 0, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k1_2p= ";  rk = 0.120442708333329_DP;  !// calulated from N1001; 0.12044270784428178 - fromN221
  res = yk_lcr_rk(f2p, f2p, f2p, f2p, 1, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

  help = "R_k2_2p= ";  rk = 0.08789062499999964_DP;  !// calulated from N1001; 0.08789062488574255 - fromN221
  res = yk_lcr_rk(f2p, f2p, f2p, f2p, 2, lcr);   call dbg_dp(TRIM(help), res);
  call dbg_assert_abs_LE(TRIM(help), res - rk, rkErr);

!  // 2p, k=2
  call atom_y2v_2p(tk, lcr%rv);                call vec_dbg("tY_k2_2p=", tk);
  call yk_lcr(ck, f2p, f2p, 2, lcr);           call vec_dbg("cY_k2_2p=", ck);
  call vec_assert_abs_diff_LE("Y_k2_2p = ", ck, tk, yzErr); 

  call mem_del_vec(f, f2);
  call mem_del_vec(f2p);
end SUBROUTINE
!___________________________________________________________
end module
!___________________________________________________________
!
!  VERSION-2
! java code YkLcrV2
!============================================================== 
!package atom.wf.lcr.yk;
!import atom.wf.lcr.TransLcrToR;
!import atom.wf.lcr.WFQuadrLcr;
!import math.Mathx;
!import math.func.FuncVec;
!import math.vec.Vec;
!import math.vec.grid.StepGrid;
!
!import javax.utilx.log.Log;
!/**
! * Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com 23/07/12, 11:47 AM
! */
!public class YkLcrV2 {
!public static Log log = Log.getLog(YkLcrV2.class);
!private final Vec wf;
!private final Vec wf2;
!private final int k;
!private final WFQuadrLcr quadr;
!private final double[] cr2;//(c+r)^2
!private final double[] cr;//(c+r)
!private final double[] r;
!private final double[] divR;
!private final double[] yDivR;
!private final Vec vR;
!
!private final FuncVec fvRes; // tmp working space AND results
!private final double[] res; // tmp working space AND results
!private final double h;
!
!
!public YkLcrV2(final WFQuadrLcr quadr, final Vec wf, final Vec wf2, final int k) {
!  this.k = k;
!  this.wf = wf;
!  this.wf2 = wf2;
!  this.quadr = quadr;
!  vR = quadr.getR();
!  TransLcrToR xToR = quadr.getLcrToR();
!  this.r = xToR.getArr();
!  this.divR = xToR.getDivR().getArr();
!  this.yDivR = xToR.getCRDivR().getArr();
!  this.cr2 = xToR.getCR2().getArr();
!  this.cr = xToR.getCR().getArr();
!
!  h = ((StepGrid) xToR.getX()).getGridStep();
!  fvRes = new FuncVec(vR);  // working space
!  res = fvRes.getArr();
!}
!public FuncVec calcZk() {
!  return calcZkV2();
!}
!public FuncVec calcYk() {
!  calcZkV2(); // result is in fvRes;
!  return calcYkV2();
!}
!
!private FuncVec calcYkV2() {
!  if (r[0] != 0) {
!    throw new IllegalArgumentException(log.error("r[0] != 0"));
!  }
!  int K1 = k + 1;
!  double C = 2 * k + 1;
!  double h2 = h / 2.0;
!  double h3 = h / 3.0;
!  double h90 = h / 90. * C;// NOTE *C
!  double h114 = 114.0 * h90;
!  double a34 = 34. * h90;
!
!  int last = res.length - 1;  // starting from last
!  double f1 = yDivR[last]   * res[last];
!  double f2 = yDivR[last-1] * res[last-1];
!  double f3 = yDivR[last-2] * res[last-2];
!  double f4 = yDivR[last-3] * res[last-3];
!  for (int i = last - 2; i >= 2; i--) {
!    double f5 = yDivR[i - 2] * res[i - 2];
!    double rM1 = r[i - 1]; // NOTE K1 not K below
!    double a3 = Mathx.pow(rM1 * divR[i + 2], K1) * h90;
!    double a2 = Mathx.pow(rM1 * divR[i + 1], K1);
!    double a1 = Mathx.pow(rM1 * divR[i + 0], K1) * h114;
!    double ai = Mathx.pow(rM1 * divR[i - 2], K1) * h90;
!    //      YK(M) = YK(M+2)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)        AATK4078
!    res[i - 1] = res[i + 1] * a2 + (a1 * f3 + a34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
!    f1 = f2;
!    f2 = f3;
!    f3 = f4;
!    f4 = f5;
!  }
!  res[0] = 0;
!  return fvRes;
!}
!
!
!private FuncVec calcZkV2() {       //log.setDbg();
!  double h2 = h / 2.0;
!  double h3 = h / 3.0;
!  double h90 = h / 90.;
!  double h114 = 114.0 * h90;
!  double a34 = 34. * h90;
!
!  res[0] = 0;
!  int n = res.length;
!
!  double f1 = cr2[0] * wf.get(0) * wf2.get(0);
!  if (f1 != 0) {
!    throw new IllegalArgumentException(log.error("(f1 != 0)"));
!  }
!  double f2 = cr2[1] * wf.get(1) * wf2.get(1);
!  double f3 = cr2[2] * wf.get(2) * wf2.get(2);
!  double f4 = cr2[3] * wf.get(3) * wf2.get(3);
!  double a, b, intgl;
!
!  // NOTE!!! Integrating by r (not X)
!  // y(r) = a r^b + c
!  // f1 = c
!  // f2 = a r2^b + c
!  // f3 = a r3^b + c
!  // f2 - f1 = a r1^b
!  // f3 - f1 = a r2^b
!  // (f3-f1)/(f2-f1)= (r3/r2)^b;
!  double r2 = r[1];
!  double r3 = r[2];
!  double rk2 = Math.pow(r2, k);
!  double rk3 = Math.pow(r3, k);
!  double f21 = f2 * rk2 / cr[1];
!  double f31 = f3 * rk3 / cr[2];
!  b = Math.log(f31 / f21) / Math.log(r3 / r2);
!  if (b <= 0) {
!    throw new IllegalArgumentException(log.error("(b <= 0)"));
!  }
!  double hb = Math.pow(r2, b);
!  a = f21 / hb;
!  intgl = a *hb / (b+1) * r2;
!  res[1] = intgl / rk2;
!
!  hb = Math.pow(r3, b);
!  intgl = a *hb / (b+1) * r3;
!  res[2] = intgl / rk3;
!
!  // START DEBUG CODE ====================================
!  //
!  // Simpson rule
!  // y(x) = a x^2 + b x + c
!  // f1 = c
!  // f2 = a h^2 + b h + c
!  // f3 = 4 a h^2 + 2 b h + c
!  //  double a = (f3 - 2. * f2 - c) / (2. * h * h);
!  //  double b = (4. * f2 - f3 + 3c) / (2. * h);
!  // int_0^h y(x) dx = a / 3. * h^3 + b h^2 / 2. + c h;
!  //  a = (f3 - 2. * f2 + f1) / 2.;
!  //  b = (4. * f2 - f3 - 3. * f1) / 2.;
!  //  intgl = (a / 3.  + b / 2. + f1) * h;
!  //  res[1] = res[0] + intgl;
!  //  res[1] = h2 * f2;   // trapezoidal
!  //
!  //  // Simpson rule
!  //  a = Mathx.pow(r[1] * divR[2], k);
!  //  //YK(3) = YK(1)*A2 + H3*(F3 + D4*A*F2 + A2*F1)
!  //  res[2] = h3*(f3 + 4.*a*f2);// NOTE!!! r[0]=Zk[0]=a2=f1=0 in the LCR grid
!  //
!  // END DEBUG CODE ====================================
!
!  for (int i = 4; i < n; i++) {
!    double f5 = cr2[i] * wf.get(i) * wf2.get(i);//      F5 = (RR(M)*P(M,I))*P(M,J)                                        AATK4125
!
!    // Zim = (ri/rim)^k Zi + I^im_i FF y^2 * (r/rim)^k dx;  where ri=r_i, rim=r_(i+m)
!    // the same as in Froese-Fischer but replacing
!    //  exp(-mhk) with (ri/rim)^k;   and
!    //  r*exp(k(x-xim)) with  y^2*(r/rim)^k
!
!    // Replace exp(-mhk) with (ri/rim)^k;  where m=2
!    //      EH = DEXP(-h)   // from SUBROUTINE INIT
!    //      A = EH**K
!    double ovR = divR[i - 1];
!    double a3 = Mathx.pow(r[i - 4] * ovR, k) * h90; //      A3 = A2*A*H90                                                     AATK4112
!    double a2 = Mathx.pow(r[i - 3] * ovR, k);       //      A2 = A*A                                                          AATK4110
!    double a1 = Mathx.pow(r[i - 2] * ovR, k) * h114;//      AN = 114.D0*A*H90                                                 AATK4114
!    double ai = Mathx.pow(r[i - 0] * ovR, k) * h90; //      AI = H90/A                                                        AATK4113
!
!    // From Froese Fischer etal Computational atomic structure 2000
!    // DELTA^2 F = F2 - 2*F3 + F4
!    // DELTA^4 F = F1 - 4*F2 + 6*F3 - 4*F4 + F5
!    // INTL F(x)dx = h(F3 + DELTA^2 F/3 - DELTA^4 F/90) + O(h^7) =
!    //  = h/90 [114*F3 + 34*(F2+F4) - F1 - F5]
!    //
!    //      YK(M-1) = YK(M-3)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)
!    res[i-1] = res[i-3] * a2 + (a1 * f3 + a34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
!    // NOTE: missing i=n-1
!    f1 = f2;
!    f2 = f3;
!    f3 = f4;
!    f4 = f5;
!  }
!  // NOTE: missing i=n-1
!  double norm = 0;
!  double rkn2 = 1;
!  double rkn3 = 1;
!  if (k == 0) {  // correct from last two points
!    norm = quadr.calcInt(wf, wf2);
!  }
!  else {
!    norm = quadr.calcInt(wf, wf2, quadr.getPowR(k));
!    rkn2 = 1. / quadr.getPowR(k).get(n - 2);
!    rkn3 = 1. / quadr.getPowR(k).get(n - 3);
!  }
!  double corr = norm * rkn2 - res[n - 2];      //log.info("corr=", corr);
!  for (int i = n-2; i > 0; i -= 2) { // do not correct i=0
!    res[i] += corr;
!  }
!  corr = norm * rkn3 - res[n - 3];             //log.info("corr=", corr);
!  for (int i = n-3; i > 0; i -= 2) { // do not correct i=0
!    res[i] += corr;
!  }
!
!  //last point
!  double a1 = Mathx.pow(r[n-2] * divR[n-1], k);
!  res[n-1] = a1 * res[n-2];
!  return fvRes;
!}
!}
!
!
!
!  VERSION-1 java conversion 
!=================================================
!package atom.wf.lcr.yk;
!import atom.wf.lcr.TransLcrToR;
!import atom.wf.lcr.WFQuadrLcr;
!import math.func.intrg.IntgInftyPts7;
!import math.func.intrg.IntgPts7;
!import math.vec.Vec;
!import math.vec.VecDbgView;
!import math.vec.grid.StepGrid;
!import math.func.FuncVec;
!import math.Mathx;
!import math.interpol.PolynIntrp;
!
!import javax.utilx.log.Log;
!// YkLogR(r)=r I^oo_0 r<^k/r>^(k+1) s^2 ds Rn(s)Rn'(s)
!// R(r)=sqrt(y)/r * P(r); where y(r)=c+r
!// YkLogR(r) = Zk(r) + I^oo_r (r/s)^(k+1) yds Pn(s) Pn'(s)
!// Zk(r) = I^r_0 (s/r)^k yds Pn(s) Pn'(s)
!//
!// dZk/dr = -k/r Zk + PPy
!// dYk/dr = (k+1)/r YkLogR - (2k+1)/r Zk
!// Zk(0)=0; YkLogR(r->oo)->Zk(r)
!// dZ/dr=dZ/dx 1/y; since x=ln(y)
!//
!// dZ/dx = -ky/r atomZ + PPy^2
!// dY/dx = (k+1)y/r Y - (2k+1)y/r atomZ
!//
!// integrating factor: r^k    for atomZ
!// d(r^k atomZ)/dx = k r^(k-1) dr/dx atomZ + r^k dZ/dx;  since dr/dx=y
!// d(r^k atomZ)/dx = FF y^2 r^k
!//
!// Zim = (ri/rim)^k Zi + I^im_i FF y^2 * (r/rim)^k dx;  where ri=r_i, rim=r_(i+m)
!// the same as in Froese-Fischer but replacing
!//  exp(-mhk) with (ri/rim)^k;   and
!//  r*exp(k(x-xim)) with  y^2*(r/rim)^k
!//
!// integrating factor: r^-(k+1)    for Y
!// d(r^-(k+1) Y)/dx = -(k+1) r^(k+1) y/r atomZ + r^-(k+1) dZ/dx;  since dr/dx=y
!// d(r^-(k+1) Y)/dx = -(2k+1)y/r r^-(k+1) atomZ
!//
!// Yim = (rim/ri)^(k+1) Yi + I^im_i atomZ y/r (rim/r)^(k+1) dx;  where ri=r_i, rim=r_(i+m)
!// the same as in Froese-Fischer but replacing
!//  exp(-mh(k+1)) with (ri/rim)^(k+1);   and
!//  exp(-(k+1)*(x-xi)) with  y/r*(r/rim)^(k+1)
!//
!public class YkLcrV1 {
!public static Log log = Log.getLog(YkLcrV1.class);
!private final static double[] yTmp = new double[10]; // tmp array for zk_OLD
!private final Vec wf;
!private final Vec wf2;
!private final double[] cr2;//(c+r)^2
!private final double[] cr;//(c+r)
!private final double[] r;
!private final Vec vR;
!private final double[] divR;
!private final double[] yDivR;
!private final int k;
!private final int MX;
!private final double H;
!private final double H2;
!private final double H3;
!private final TransLcrToR xToR;
!private final WFQuadrLcr quadr;
!
!// NOTE! THIS IS optimization. Stop re-calculating the same values
!private static FuncVec zF;  // z-function inside int_0^r
!private static FuncVec rK;    // 1/r^k
!private static FuncVec yF;  // y-function inside int_0^r
!private static FuncVec rK1;    // r^(k+1)
!private static FuncVec last_f;
!private static FuncVec last_f2;
!private static Vec last_r;
!private static int last_K;
!private boolean useLast = false;
!
!public YkLcrV1(final WFQuadrLcr quadr, final Vec wf, final Vec wf2, final int K) {
!  this.k = K;
!  this.wf = wf;
!  this.wf2 = wf2;
!  this.quadr = quadr;
!  xToR = quadr.getLcrToR();
!  this.r = xToR.getArr();
!  this.vR = xToR;
!  this.cr2 = xToR.getCR2().getArr();
!  this.cr = xToR.getCR().getArr();
!  this.divR = xToR.getDivR().getArr();
!  this.yDivR = xToR.getCRDivR().getArr();
!  if (!(xToR.getX() instanceof StepGrid)) {
!    String mssg = "YkLogR can only work with StepGrid";
!    throw new IllegalArgumentException(log.error(mssg));
!  }
!  this.H = ((StepGrid) xToR.getX()).getGridStep();
!  this.H2 = H / 2.0;
!  this.H3 = H / 3.0;
!  this.MX = (Math.min(wf.size(), wf2.size()) / 2) * 2;//      MX = (MIN0(MAX(I),MAX(J))/2)*2                                    AATK4123
!  // MX is used to trace simpson's rule errors
!  //EH = DEXP(-H)   // from SUBROUTINE INIT
!
!  useLast = checkUseLast();
!}
!private boolean checkUseLast() {
!  if (last_f == wf && last_f2 == wf2 &&  last_K == k &&  last_r  == vR)
!    return true;
!  return false;
!}
!public FuncVec calcZk() {
!//  return calcZk_DBG();
!  return calcZk_v1();
!}
!public FuncVec calcYk() {
!//  return calcYk_DBG();
!  return calcYk_v1();
!}
!public FuncVec calcYk_DBG() {
!  FuncVec zk = calcZk_DBG();
!//  FuncVec zk = calcZk_v1();
!  return calcYk_DBG(zk);
!}
!public FuncVec calcYk_v1() {
!  FuncVec zk = calcZk_v1();
!//  FuncVec zk = calcZk_DBG();
!  return calcYk_v1(zk);
!}
!public FuncVec calcYk_DBG(FuncVec zk) {   //log.setDbg();
!  loadYFuncs(zk);
!  // integrate backwards
!  FuncVec res = new IntgInftyPts7(yF);    log.info("IntgInftyPts7(zF)=", new VecDbgView(res));
!
!  // boundary Y_k(r-->oo) = Z_k(r)
!  double corr = zk.getLast() / rK1.getLast();
!  res.add(corr - res.getLast());
!
!  res.multSelf(rK1);
!  return res;
!}
!public FuncVec calcZk_DBG() {   //log.setDbg();
!  loadZFuncs();
!
!  // IntgInftyPts7 is not good for ln(r)-type grid!!!!!!!!
!  // delta_r is much smaller for small r's than for large r's,
!  // so it is mor accurate to integrate from start
!//  FuncVec res = new IntgInftyPts7(zF);    log.info("IntgInftyPts7(zF)=", new VecDbgView(res));
!//  res.add(-res.getFirst());
!
!  //NOTE!!  HERE  IntgPts7 is much better than calcFuncInt_DEV
!  FuncVec res = new IntgPts7(zF);    log.info("IntgPts7(zF)=", new VecDbgView(res));
!//  FuncVec res = quadr.calcFuncInt_DEV(zF);    log.info("IntgPts7(zF)=", new VecDbgView(res));
!
!  res.multSelf(rK);
!  return res;
!}
!private void loadZFuncs() {
!  zF = new FuncVec(xToR.getX());
!  double[] z = zF.getArr();
!  rK = new FuncVec(vR);
!  double[] rk = rK.getArr();
!  z[0] = 0;
!  rk[0] = 0;
!  for (int i = 1; i < z.length; i++) {
!    z[i] = cr2[i] * wf.get(i) * wf2.get(i);
!    double rki = Mathx.pow(r[i], k);
!    z[i] *= rki;
!    rk[i] = 1. / rki;
!  }
!  log.info("zF=", new VecDbgView(zF));
!  log.info("rK=", new VecDbgView(rK));
!}
!private void loadYFuncs(FuncVec zkFunc) {
!  yF = new FuncVec(xToR.getX());
!  double[] y = yF.getArr();
!  rK1 = new FuncVec(vR);
!  double[] rk = rK1.getArr();
!  double[] zk = zkFunc.getArr();
!  y[0] = 0;
!  rk[0] = 0;
!  double kk = -(2. * k + 1.);
!  int K1 = k + 1;
!  for (int i = 1; i < y.length; i++) {
!    y[i] = kk * (cr[i] / r[i])* zk[i];
!    double rki = Mathx.pow(r[i], K1);
!    y[i] /= rki;
!    rk[i] = rki;
!  }
!  log.info("yF=", new VecDbgView(yF));
!  log.info("rK1=", new VecDbgView(rK1));
!}
!public FuncVec calcZk_v1() {
!  FuncVec res = new FuncVec(vR);
!  double[] YK = res.getArr();
!  //      DEN = L(I) + L(J) + 3+ K                                          AATK4107
!  //      FACT = (D1/(L(I)+1) + D1/(L(J)+1))/(DEN + D1)                     AATK4108
!  //      A = EH**K                                                         AATK4109
!  //      A2 = A*A                                                          AATK4110
!  double H90 = H / 90.0;//      H90 = H/90.D0                                                     AATK4111
!  //      A3 = A2*A*H90                                                     AATK4112
!  //      AI = H90/A                                                        AATK4113
!  //      AN = 114.D0*A*H90                                                 AATK4114
!  double A34 = 34.0 * H90;//      A34 = 34.D0*H90                                                   AATK4115
!  int M = 0;
!  double f1 = cr2[M] * wf.get(M) * wf2.get(M);
!  M++;//      F1 = RR(1)*P(1,I)*P(1,J)                                          AATK4116
!  double f2 = cr2[M] * wf.get(M) * wf2.get(M);
!  M++;//      F2 = RR(2)*P(2,I)*P(2,J)                                          AATK4117
!  double f3 = cr2[M] * wf.get(M) * wf2.get(M);
!  M++;//      F3 = RR(3)*P(3,I)*P(3,J)                                          AATK4118
!  double f4 = cr2[M] * wf.get(M) * wf2.get(M);
!  M++;//      F4 = RR(4)*P(4,I)*P(4,J)                                          AATK4119
!  //YK(1) = F1*(D1 + atomZ*R(1)*FACT)/DEN                                 AATK4120
!  //YK(2) = F2*(D1 + atomZ*R(2)*FACT)/DEN                                 AATK4121
!  //YK(3) = YK(1)*A2 + H3*(F3 + D4*A*F2 + A2*F1)                      AATK4122
!  YK[0] = 0;
!  YK[1] = approxFirstZ(1, f2, f3);
!  YK[2] = approxFirstZ(2, f2, f3); // simpson does not work as good as this
!//  YK[1] = 0;
!//  YK[2] = 0;
!
!  for (M = 5; M <= MX; M++) {//      DO 8 M = 5,MX                                                     AATK4124
!    double f5 = cr2[M - 1] * wf.get(M - 1) * wf2.get(M - 1);//      F5 = (RR(M)*P(M,I))*P(M,J)                                        AATK4125
!
!    // Zim = (ri/rim)^k Zi + I^im_i FF y^2 * (r/rim)^k dx;  where ri=r_i, rim=r_(i+m)
!    // the same as in Froese-Fischer but replacing
!    //  exp(-mhk) with (ri/rim)^k;   and
!    //  r*exp(k(x-xim)) with  y^2*(r/rim)^k
!
!    // Replace exp(-mhk) with (ri/rim)^k;  where m=2
!    //      EH = DEXP(-H)   // from SUBROUTINE INIT
!    //      A = EH**K                                                         AATK4109
!    double ovR = divR[M - 2];
!    double a3 = Mathx.pow(r[M - 5] * ovR, k) * H90;//      A3 = A2*A*H90                                                     AATK4112
!    double a2 = Mathx.pow(r[M - 4] * ovR, k);//      A2 = A*A                                                          AATK4110
!    double a1 = Mathx.pow(r[M - 3] * ovR, k) * 114.0 * H90;//      AN = 114.D0*A*H90                                                 AATK4114
!    double ai = Mathx.pow(r[M - 1] * ovR, k) * H90;//      AI = H90/A                                                        AATK4113
!
!    // From Froese Fischer etal Computational atomic structure 2000
!    // DELTA^2 F = F2 - 2*F3 + F4
!    // DELTA^4 F = F1 - 4*F2 + 6*F3 - 4*F4 + F5
!    // INTL F(x)dx = h(F3 + DELTA^2 F/3 - DELTA^4 F/90) + O(h^7) =
!    //  = h/90 [114*F3 + 34*(F2+F4) - F1 - F5]
!    //
!    //      YK(M-1) = YK(M-3)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)      AATK4126
!    double YKM_3 = YK[M - 4];
!    double YKM_1 = YKM_3 * a2 + (a1 * f3 + A34 * (f4 + a2 * f2) - f5 * ai - f1 * a3);
!    YK[M - 2] = YKM_1;
!    f1 = f2;//      F1 = F2                                                           AATK4127
!    f2 = f3;//      F2 = F3                                                           AATK4128
!    f3 = f4;//      F3 = F4                                                           AATK4129
!    f4 = f5;//8     F4 = F5                                                           AATK4130
!  }
!  int M1 = MX - 1;//      M1 = MX - 1                                                       AATK4131
!//  if (K == 0) {//      IF (IABS(I-J)  +  IABS(K) .NE. 0) GO TO 2                         AATK4132
!//    //*  *****  FOR Y0(I,I) SET THE LIMIT TO 1 AND REMOVE OSCILLATIONS        AATK4134
!//    //*  *****  INTRODUCED BY THE USE OF SIMPSON'S RULE                       AATK4135
!//    int M2 = M1 - 1;//      M2 = M1 - 1                                                       AATK4137
!//    double mid = 0.5 * (YK[M1 - 1] + YK[M2 - 1]);
!//    double C1 = mid - YK[M1 - 1];//      C1 = D1 - YK(M1)                                                  AATK4138
!//    double C2 = mid - YK[M2 - 1];//      C2 = D1 - YK(M2)                                                  AATK4139
!//    for (M = 1; M <= M1; M += 2) {//      DO 3 M = 1,M1,2                                                   AATK4140
!//      YK[M - 1] += C1;//      YK(M) = YK(M) + C1                                                AATK4141
!//      YK[M] += C2;//3     YK(M+1) = YK(M+1) + C2                                            AATK4142
!//    }
!//  }
!  int NO = r.length;
!  for (M = M1 + 1; M <= NO; M++) {//2     DO 1 M = M1+1,NO                                                  AATK4143
!    double A = Mathx.pow(r[M - 2] * divR[M - 1], k);//A = EH**K                                                         AATK4109
!    YK[M - 1] = A * YK[M - 2];//         YK(M) = A*YK(M-1)                                              AATK4144
!  }//1     CONTINUE                                                          AATK4145
!  return res;
!}
!private double approxFirstZ(int idxDest, double F2, double F3) {
!  FuncVec tmp = new FuncVec(vR, yTmp);
!  int M = 0;
!  tmp.set(M, 0);
!  M++;
!  double fx[] = {0, 0, 0};
!  double rDest = r[idxDest];
!  double A = Mathx.pow(r[M] / rDest, k);
!  fx[M] = F2 * A / cr[M];
!  tmp.set(M, fx[M]);
!  M++;
!  A = Mathx.pow(r[M] / rDest, k);
!  fx[M] = F3 * A / cr[M];
!  tmp.set(M, fx[M]);
!  double b = PolynIntrp.calcPowerSLOW(tmp, 0);
!  return rDest / (b + 1) * fx[idxDest];
!}
!public FuncVec calcYk_v1(FuncVec res) {
!  double[] YK = res.getArr();
!//  FuncVec res = calcZk(); //CALL ZK(I,J,K)                                                    AATK4062
!  //A =    EH**(K+1)                                                  AATK4063
!  int K1 = k + 1;
!  double C = 2 * k + 1;//      C = 2*K+1                                                         AATK4064
!  //      A2 = A*A                                                          AATK4065
!  double H90 = C * H3 / 30.;//      H90 = C*H3/D30                                                    AATK4066
!  //      A3 = A2*A*H90                                                     AATK4067
!  //      AI = H90/A                                                        AATK4068
!  //      AN = 114.D0*A*H90                                                 AATK4069
!  double A34 = 34. * H90;//      A34 = 34.D0*H90                                                   AATK4070
!  //      MX = (MIN0(MAX(I),MAX(J))/2)*2                                    AATK4071
!  int M = MX - 1;
!  double F1 = yDivR[M] * res.get(M);
!  M--;//      F1 = YK(MX)*EH**K                                                 AATK4072
!  double F2 = yDivR[M] * res.get(M);
!  M--;//      F2 = YK(MX)                                                       AATK4073
!  double F3 = yDivR[M] * res.get(M);
!  M--;//      F3 = YK(MX-1)                                                     AATK4074
!  double F4 = yDivR[M] * res.get(M);
!  M--;//      F4 = YK(MX-2)                                                     AATK4075
!  // NOTE!!! MX-1, MX-2, and MX-3; the last three must be set to atomZ, otherwise the algorithm does not work?!
!  for (M = MX - 3; M >= 2; M--) {//      DO 9 M = MX-2,2,-1                                                AATK4076
!    double F5 = yDivR[M - 2] * res.get(M - 2);//      F5 = YK(M-1)                                                      AATK4077
!    double rM1 = r[M - 1]; // NOTE K1 not K below
!    double A3 = Mathx.pow(rM1 * divR[M + 2], K1) * H90;//      A3 = A2*A*H90                                                     AATK4112
!    double A2 = Mathx.pow(rM1 * divR[M + 1], K1);//      A2 = A*A                                                          AATK4110
!    double AN = Mathx.pow(rM1 * divR[M + 0], K1) * 114.0 * H90;//      AN = 114.D0*A*H90                                                 AATK4114
!    double AI = Mathx.pow(rM1 * divR[M - 2], K1) * H90;//      AI = H90/A                                                        AATK4113
!
!    //      YK(M) = YK(M+2)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)        AATK4078
!    YK[M - 1] = YK[M + 1] * A2 + (AN * F3 + A34 * (F4 + A2 * F2) - F5 * AI - F1 * A3);
!    F1 = F2;//      F1 = F2                                                           AATK4079
!    F2 = F3;//      F2 = F3                                                           AATK4080
!    F3 = F4;//      F3 = F4                                                           AATK4081
!    F4 = F5;//9     F4 = F5                                                           AATK4082
!  }
!  //      YK(1) = YK(3)*A2+C*H3*(F4 + D4*A*F3 + A2*F2)                      AATK4083
!  YK[0] = 0; // always zero
!  return res;
!}
!}
!
!
!=== ORIGINAL =============================================================
!/*
!*AATKHF86.  GENERAL HARTREE-FOCK PROGRAM.  C.F. FISCHER.                 AATK0000
!*REF. IN COMP. PHYS. COMMUN. 43 (1987) 355                               AATK0000
!*     ------------------------------------------------------------------AATK0001
!*                                                                       AATK0002
!*     A GENERAL HARTREE-FOCK PROGRAM                                    AATK0003
!*                                                                       AATK0004
!*     by C. Froese Fischer                                              AATK0005
!*        Vanderbilt University                                          AATK0006
!*        Nashville, TN 37235 USA                                        AATK0007
!*                                                                       AATK0008
!*     April, 1986                                                       AATK0009
!*                                                                       AATK0010
!*     ------------------------------------------------------------------AATK0011
!*                                                                       AATK0012
!*                                                                       AATK0013
!*     All comments in the program listing assume the radial function P  AATK0014
!*     is the solution of an equation of the form                        AATK0015
!*                                                                       AATK0016
!*      P" + ( 2Z/R - Y - L(L+1)/R**2 - E)P = X + T                      AATK0017
!*                                                                       AATK0018
!*     where Y is called a potential function                            AATK0019
!*           X is called an exchange function, and                       AATK0020
!*           T includes contributions from off-diagonal energy parameter AATK0021
!*                                                                       AATK0022
!*     The program uses LOG(atomZ*R) as independent variable and             AATK0023
!*                      P/SQRT(R) as dependent variable.                 AATK0024
!*     As a result all equations must be transformed as described in     AATK0025
!*     Sec. 6-2 and 6-4 of the book - ``The Hartree-Fock Method for      AATK0026
!*     Atoms'',Wiley Interscience, 1977, by Charlotte FROESE FISCHER.    AATK0027
!*     (All references to equations and sections refer to this book.)    AATK0028
!*                                                                       AATK0029
!*     This program is written as a double precision FORTRAN77           AATK0030
!*     program.  However, on computers with a word length of 48 bits or  AATK0031
!*     more it should be run as a single precision program.  Such con-   AATK0032
!*     version is facilitated through the use of IMPLICIT type declar-   AATK0033
!*     ations and the initialization of many double precision            AATK0034
!*     constants in the INIT program.                                    AATK0035
!*                                                                       AATK0036
!*                                                                       AATK0037
!*                                                                       AATK0051
!*     ------------------------------------------------------------------AATK4046
!*                 Y K F                                                 AATK4047
!*     ------------------------------------------------------------------AATK4048
!*                                                                       AATK4049
!*               k                                                       AATK4050
!*       Stores Y (i, j; r) in the array YK                              AATK4051
!*                                                                       AATK4052
!*                                                                       AATK4053
!SUBROUTINE YKF(I,J,K,REL)                                         AATK4054
!IMPLICIT DOUBLE PRECISION(A-H,O-atomZ)                                AATK4055
!PARAMETER (NOD=220,NWFD=20)                                       AATK4056
!COMMON /PARAM/H,H1,H3,CH,EH,RHO,atomZ,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID  AATK4057
!:   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD      AATK4058
!COMMON /RADIAL/R(NOD),RR(NOD),R2(NOD),P(NOD,NWFD),YK(NOD),        AATK4059
!:   YR(NOD),X(NOD),AZ(NWFD),L(NWFD),MAX(NWFD),N(NWFD)              AATK4060
!LOGICAL REL                                                       AATK4061
!CALL ZK(I,J,K)                                                    AATK4062
!A =    EH**(K+1)                                                  AATK4063
!C = 2*K+1                                                         AATK4064
!A2 = A*A                                                          AATK4065
!H90 = C*H3/D30                                                    AATK4066
!A3 = A2*A*H90                                                     AATK4067
!AI = H90/A                                                        AATK4068
!AN = 114.D0*A*H90                                                 AATK4069
!A34 = 34.D0*H90                                                   AATK4070
!MX = (MIN0(MAX(I),MAX(J))/2)*2                                    AATK4071
!F1 = YK(MX)*EH**K                                                 AATK4072
!F2 = YK(MX)                                                       AATK4073
!F3 = YK(MX-1)                                                     AATK4074
!F4 = YK(MX-2)                                                     AATK4075
!DO 9 M = MX-2,2,-1                                                AATK4076
!F5 = YK(M-1)                                                      AATK4077
!YK(M) = YK(M+2)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)        AATK4078
!F1 = F2                                                           AATK4079
!F2 = F3                                                           AATK4080
!F3 = F4                                                           AATK4081
!9     F4 = F5                                                           AATK4082
!YK(1) = YK(3)*A2+C*H3*(F4 + D4*A*F3 + A2*F2)                      AATK4083
!IF (.NOT.REL) RETURN                                              AATK4084
!C = C*FINE                                                        AATK4085
!DO 10 M = 1,MX                                                    AATK4086
!YK(M) = YK(M) + C*P(M,I)*P(M,J)                                   AATK4087
!10    CONTINUE                                                          AATK4088
!RETURN                                                            AATK4089
!END                                                               AATK4090
!*                                                                       AATK4091
!*AATKHF86.  GENERAL HARTREE-FOCK PROGRAM.  C.F. FISCHER.                 AATK0000
!*REF. IN COMP. PHYS. COMMUN. 43 (1987) 355                               AATK0000
!*     ------------------------------------------------------------------AATK0001
!*                                                                       AATK0002
!*     A GENERAL HARTREE-FOCK PROGRAM                                    AATK0003
!*                                                                       AATK0004
!*     by C. Froese Fischer                                              AATK0005
!*        Vanderbilt University                                          AATK0006
!*        Nashville, TN 37235 USA                                        AATK0007
!*                                                                       AATK0008
!*     April, 1986                                                       AATK0009
!*                                                                       AATK0010
!*     ------------------------------------------------------------------AATK0011
!*                                                                       AATK0012
!*                                                                       AATK0013
!*     All comments in the program listing assume the radial function P  AATK0014
!*     is the solution of an equation of the form                        AATK0015
!*                                                                       AATK0016
!*      P" + ( 2Z/R - Y - L(L+1)/R**2 - E)P = X + T                      AATK0017
!*                                                                       AATK0018
!*     where Y is called a potential function                            AATK0019
!*           X is called an exchange function, and                       AATK0020
!*           T includes contributions from off-diagonal energy parameter AATK0021
!*                                                                       AATK0022
!*     The program uses LOG(atomZ*R) as independent variable and             AATK0023
!*                      P/SQRT(R) as dependent variable.                 AATK0024
!*     As a result all equations must be transformed as described in     AATK0025
!*     Sec. 6-2 and 6-4 of the book - ``The Hartree-Fock Method for      AATK0026
!*     Atoms'',Wiley Interscience, 1977, by Charlotte FROESE FISCHER.    AATK0027
!*     (All references to equations and sections refer to this book.)    AATK0028
!*                                                                       AATK0029
!*     This program is written as a double precision FORTRAN77           AATK0030
!*     program.  However, on computers with a word length of 48 bits or  AATK0031
!*     more it should be run as a single precision program.  Such con-   AATK0032
!*     version is facilitated through the use of IMPLICIT type declar-   AATK0033
!*     ations and the initialization of many double precision            AATK0034
!*     constants in the INIT program.                                    AATK0035
!*                                                                       AATK0036
!*                                                                       AATK0037
!*     ------------------------------------------------------------------AATK4092
!*                 atomZ K                                                   AATK4093
!*     ------------------------------------------------------------------AATK4094
!*                                                                       AATK4095
!*               k                                                       AATK4096
!*       Stores atomZ (i, j; r) in the array YK.                             AATK4097
!*                                                                       AATK4098
!*                                                                       AATK4099
!SUBROUTINE ZK(I,J,K)                                              AATK4100
!PARAMETER (NOD=220,NWFD=20)                                       AATK4101
!IMPLICIT DOUBLE PRECISION(A-H,O-atomZ)                                AATK4102
!COMMON /PARAM/H,H1,H3,CH,EH,RHO,atomZ,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID  AATK4103
!:   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD      AATK4104
!COMMON /RADIAL/R(NOD),RR(NOD),R2(NOD),P(NOD,NWFD),YK(NOD),        AATK4105
!:   YR(NOD),X(NOD),AZ(NWFD),L(NWFD),MAX(NWFD),N(NWFD)              AATK4106
!DEN = L(I) + L(J) + 3+ K                                          AATK4107
!FACT = (D1/(L(I)+1) + D1/(L(J)+1))/(DEN + D1)                     AATK4108
!A = EH**K                                                         AATK4109
!A2 = A*A                                                          AATK4110
!H90 = H/90.D0                                                     AATK4111
!A3 = A2*A*H90                                                     AATK4112
!AI = H90/A                                                        AATK4113
!AN = 114.D0*A*H90                                                 AATK4114
!A34 = 34.D0*H90                                                   AATK4115
!F1 = RR(1)*P(1,I)*P(1,J)                                          AATK4116
!F2 = RR(2)*P(2,I)*P(2,J)                                          AATK4117
!F3 = RR(3)*P(3,I)*P(3,J)                                          AATK4118
!F4 = RR(4)*P(4,I)*P(4,J)                                          AATK4119
!YK(1) = F1*(D1 + atomZ*R(1)*FACT)/DEN                                 AATK4120
!YK(2) = F2*(D1 + atomZ*R(2)*FACT)/DEN                                 AATK4121
!YK(3) = YK(1)*A2 + H3*(F3 + D4*A*F2 + A2*F1)                      AATK4122
!MX = (MIN0(MAX(I),MAX(J))/2)*2                                    AATK4123
!DO 8 M = 5,MX                                                     AATK4124
!F5 = (RR(M)*P(M,I))*P(M,J)                                        AATK4125
!YK(M-1) = YK(M-3)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)      AATK4126
!F1 = F2                                                           AATK4127
!F2 = F3                                                           AATK4128
!F3 = F4                                                           AATK4129
!8     F4 = F5                                                           AATK4130
!M1 = MX - 1                                                       AATK4131
!IF (IABS(I-J)  +  IABS(K) .NE. 0) GO TO 2                         AATK4132
!*                                                                       AATK4133
!*  *****  FOR Y0(I,I) SET THE LIMIT TO 1 AND REMOVE OSCILLATIONS        AATK4134
!*  *****  INTRODUCED BY THE USE OF SIMPSON'S RULE                       AATK4135
!*                                                                       AATK4136
!M2 = M1 - 1                                                       AATK4137
!C1 = D1 - YK(M1)                                                  AATK4138
!C2 = D1 - YK(M2)                                                  AATK4139
!DO 3 M = 1,M1,2                                                   AATK4140
!YK(M) = YK(M) + C1                                                AATK4141
!3     YK(M+1) = YK(M+1) + C2                                            AATK4142
!2     DO 1 M = M1+1,NO                                                  AATK4143
!YK(M) = A*YK(M-1)                                              AATK4144
!1     CONTINUE                                                          AATK4145
!END                                                               AATK4146
!*/
